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Abstract 

The reliability of surface-based electrical resistivity tomography (ERT) for quan- 
tifying resistivities for shallow subsurface water processes is analysed. A method 
comprising numerical simulations of water movement in soil and forward-inverse 
modeling of ERT surveys for two synthetic data sets is presented. Resistivity con- 
trast, e.g. by changing water content, is shown to have large influence on the resis- 
tivity quantification. 

An ensemble and clustering approach is introduced in which ensembles of 50 
different inversion models for one data set are created by randomly varying the 
parameters for a regularisation based inversion routine. The ensemble members are 
sorted into five clusters of similar models and the mean model for each cluster 
is computed. Distinguishing persisting features in the mean models from singular 
artifacts in individual tomograms can improve the interpretation of inversion results. 

Especially in presence of large resistivity contrasts in high sensitivity areas, the 
quantification of resistivities can be unreliable. The ensemble approach shows that 
this is an inherent problem present for all models inverted with the regularisation 
based routine. The results also suggest that the combination of hydrological and 
electrical modeling might lead to better results. 
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1 1 Introduction 



The quantification of water content by geophysical methods is an important 
focus of hydrogeophysical research. Surface based electrical resistivity tomog- 
raphy (ERT) is a promising method, because it is non-intrusive and can cover 
large surface areas quickly, while it might also be permanently installed for au- 
tomated monitoring purposes. The development of inversion software for the 



processing of measured (apparent) resistivit ies to mode 



s of tru e resistivity has 



20041 ). Consequently, 



made fast and extensive surveys possible (iDaily et al.l . 
assessing the reliability of ERT for quantifying soil water content is a currently 
active research field. 

ERT has successfully been used in a num ber of different applications, e.g . 



12 in bo rehole surveys of tracer experi ments ( 



2OO2I ) or in laboratory experiments (IBinley et al 



Slater 



et al 



1996 



2000 


■ 


Kemna et al. 


Slater et al. 




2002) 



It has also been apphed in surface-based surveys of the vadose zone (e.g 



Daily and Ramirez 



200l|). 



19921 ) and of groundwater flow after heavy rain (jSuzuki and Higashi 



Because choice of measurement configuration and inversion parameters may 
have significant influence on the survey re sults, improving the qua lity of ERT 



19 surveys has been an intense research topic. 



Dahlin and Zhoul (12004 ) have com- 
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pared 10 different electrod e arrays for 2D 



21 using syntlietic data sets. 



Stummer et al 



survey s and assessed tlieir quality 



( I2OO4I ) have developed algorithms 



22 to calculate optimal electrode arrays that provide as much i nformation on the 



23 subsurface as possible. 



2003 



Oldenborger et al. 



Sioedahl et al 



he effects of measurement errors ( 



2OO5I ) and geometry (ILoke 



20061 ) and inversion parameters (ICarle et al. 



2000 



Zhou and Dahlin 



Hennig et al.l. 12005 



1999 



Rings et al. 



20081 ) on the surveys have been studied. 

Geophysical methods cannot directly determine hydrological properties like 
soil water content. They must be deducted using a general or calibrated rela- 
tionship between the attribute of interest and the property available through 
geophysical measurements. In the case of ERT, the resistivities of the subsur- 
face are related to wate r cont ent by a generic petrophysical relation; usually 
the equation by Archie ( 1942! ). The resistivities, again, are not readily avail- 



able from surface-based ERT surveys, but must be obtained from the measured 
apparent resistivities via inversion. The most widespread inversion methods 
rely on regularised least-squares minimisation to find the smoothest model 
of resistivities that gives a model response closest to the measured apparent 
resistivities. 

Even assuming that the petrophysical relation between resistivity and water 
content is known, the resistivity models are non-unique and have likely been 
affected by the inversion process. The sensitivity of tomographic surveys plays 
a major role in the retrieval of subsurface characteristics, e.g. for surface-based 
ERT the sensitivity decreases with depth. Low sensiti vity areas (but no t only 



Rings et al. 



20081 ). The 



those) can often be plagued by inversion artifacts (e.g. 
inversion process and the choice of inversion parameters, e.g. the regularisa- 
tion parameters, determine how well the inverted model will reproduce the 



46 real distribution. However, some of the parameter choices can not reliably be 
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47 based upon obse r vation , but must be fitted or depend on experience. 



Day-Lewis et al. 



(12004 



20051 ) refer to the loss of information caused by the 

49 inversion process, lack of sufficient prior information and survey geometry as 

50 'correlation loss'. They developed a method to compute the correlation loss as 

51 a function of the influencing factors. This allows an analytical integration of 

52 these factors into geostatistical analyses of quantitative hydrological field sur 



53 veys, but needs a priori knowledge of covariance models. 



Singha and Gorelick 



54 ( I2OO6I ) suggest a nonstationary estimation approach that uses numerical simu- 

55 lations of transport and electrical current fiow to deduct apparent petrophysi- 

56 cal relations. These methods modify the translation from the inverted models 

57 by adjusting the petrophysical relation but require either a priori knowledge 

58 or are computationally intensive. 

59 To assess the quahty of ERT-based water content quantification, the complete 

60 processing chain including the inversion process, the petrophysical relation 

61 and numerical simulations of the soil water movement has to be evaluated. 

62 This study introduces a combined approach using soil hydraulic simulations 

63 and ensemble building of inverted models to estimate the uncertainty inherent 

64 in typical applications of ERT for water content quantification. 



65 2 Methods 



66 To evaluate the inversion process, a forward-inverse cycle approach is used. In 

67 numerous applications and studies, forward modehng of synthetic data sets 

68 has been used to gain additi onal insight and confidence into measurements 



69 


and the inversion process (e.g. 


Lo 


ce and Dahlin. 


2OO2I: 


Godio and Naldi. 


2003 


70 


Hauck and Yonder Muehll. 


2003; 


Loke et al.. 


2OO3I: 


Nffuven et al.. 


2005, 


2007 



Rings et all l2005l ). Forward modeling routines are applied to synthetic data 

72 sets obtained from simulations of soil water movement. For two cases studies, 

73 the approach is used to discuss how slight variations in the soil structure 

74 influence the resistivity retrieval, and thereby the water content retrieval. 

75 The second part of the study proposes an ensemble approach which allows an 

76 overview of the possible range of inverted models, improves the analysis and 

77 enables general assertions about how well a given model can be characterised 

78 through the chosen inversion process. 

79 In the following, each methodological step of the methods will be shortly 

80 introduced, further discussion will illustrate how these steps can be applied to 

81 create and analyse two synthetic data sets. 

82 The forward-inverse cycle consists of three steps: 



83 (1) Simulation of water movement in soil: A model with specific soil structure 

84 is generated for numerical simulation of water movement. The movement 

85 of a water front, caused by infiltrating rainfall, is simulated over time. 

86 Characteristic states of water percolation are identified (starting with a 

87 completely dry soil) and a simplified distribution of water content for 

88 each state is extracted. 

89 (2) Generic resistivity model: A generic resistivity model mirroring the soil 

90 structure from (1) is created. 

91 • For a model representing a dry state (no water content), resistivities are 

92 assigned based on typical values known from laboratory measurements 

93 and/or literature. 

94 • For states of water percolation, changes in water content can be calcu- 

95 lated using the water content distribution from (1). They can be trans- 

96 ferred into resistivity changes by applying a petrophysical relation, e.g. 
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97 the equation by lArchid (jl942l ). 

98 • A finite-element based forward modeling routine transfers tlie generic 

99 resistivity models into model responses (sets of apparent resistivities) 

100 that correspond to the data that would have been recorded by field 

101 surveys. Random noise is added to simulate field measuring conditions. 

102 (3) Resistivity inversion: The apparent resistivities are inverted using a suit- 

103 able inversion scheme. The most widespread inversion schemes include 

104 smoothness constrained (L2-norm) methods and robust (Ll-norm) schemes 

105 which are preferable if sharp layer boundaries are present. The forward- 

106 inverse cycle is completed by comparing and evaluating the generic and 

107 inverted model of resistivities. 

108 The ensemble method comprises two steps: 

109 (1) Ensemble generation: For each data set, an ensemble of 50 different in- 
verted models is created by varying the inversion parameters and/or the 
inversion scheme. The parameter set is chosen randomly from a parameter 
space constrained to physically meaningful parameter sets. 

(2) Clustering: A clustering algorithm is used to group similar models of the 
ensemble. Cluster members can be averaged to simplify the analysis of 
the ensemble. 



116 2.1 Forward-inverse cycle 

111 The application of this methodology was governed by the available software 

118 codes for modeling and inversion. This section discusses how the steps were 

119 specifically realised to create and analyse two synthetic data sets. 
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120 2.1.1 Simulation of water movement in soil 



121 A numerical simulation of water movement was used to ensure that realistic 

122 distributions of water content (and thus resistivity) were used in this study. 

123 If a continuously connecte d air pha s e is a. ssumed, the equation of motion for 



124 water in soil was given by 



Richards! (Il93ll ) as: 



d_ 

dt 



126 with volumetric water content 6^, hydraulic conductivity K^j, matric potential 

127 \E'm, density of water and gravitational acceleration g. To solve Eq. [T]for 

128 water content, the material properties have to be given that connect 6'^;, 

129 and "^m- Usually, the soil- water characteristic 6'^(\E'm) and the conductivity 

130 Ku,{9w) are parameterised. 



131 The m ost widely used parameterisation for the soil- water characteristics (Ivan Genuchtenl . 



132 Il980l ). written in terms of water saturation S = (6 — 6^)/ {dg — Or) with residual 

133 water content Or, saturated volumetric water content dg and hydraulic head 

134 hm = "^m/iQwg), is 



S{hr, 



(2) 



136 with the scaling factor a, which is related to the air-entry value and the 

137 parameter u connected to the pore size distribution. The hydrauli c cond uc- 



138 tivity is characterised by applying the parameterisa tions of 



139 concise overview of the soil physics is given e.g. by 



StephensI fll996h 



140 Equa tionUwas solved numerically using the HYDRUS software (jSimunek et al. 



MualemI (119761 ). A 



20061 ). By defining time- variable precipitation and evaporation rates as atmo- 



142 spheric boundary conditions, changes in the hydraulic head hm and thus water 
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20081 ). In 



143 movement are induced. 

144 The simulations were conducted with models represen ting a two-layered soil 

145 representative of a site used in previous field studies (IRings et al.l . 

146 addition to an atmospheric boundary, a seepage boundary on the bottom al- 

147 lowed water to leave the domain. From the simulations, characteristic states 

148 of a water front infiltrating the domain were identified. Generally, beyond the 

149 dry state, characteristic states should be chosen at times when the water con- 

150 tent distribution has changed significantly, e.g. when one layer has become 

151 completely saturated. 



153 2.1.2 Generic resistivity model 



154 The transfer from water sa t uratio n values 5* to electrical resistivity p is given 



155 by the equations of 



Archid (Il942l ). Here the quotient form is applied given by 



Pi_ 



(3) 



157 where it is assumed that two measurements of the same soil at time steps i, j 

158 differing only in water saturation are connected by the saturation exponent 

159 n. n is near 2 for an organic overburden and in the range of 1.01 to 2.7 for 



160 unconsolidated sands (lUlrich and Slater 



20041 ) 



161 A generic model of resistivities was constructed calculate the response (mea- 

162 surement data) an actual ERT survey would have retrieved. We simulated 

163 field conditions by superimposing 3 % random noise on the resulting apparent 

164 resistivity data set. 
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165 2.1.3 Inversion of apparent resistivities 



166 The forward- inverse cycle is completed by inverting the simulated measure- 

167 ment data. Generic and inverted models can then be compared and the dis- 

168 crepancies analysed. 

169 A robust inversion scheme by 



Loke et al. 



(120031 ). which is usually employed 



whereever sharp layer boundaries are expected, w as chosen. It is implemented 



171 as an iteratively reweighted least-squares method (IWolke and Schwetlick 

172 in the software RES2DINV: 



19881) 



(4) 



174 Here Jj are Jacobian matrices of partial derivatives for the i-t h iteration, W is 



175 a rou ghness filter using a first-order finite-difference operator (jdeGroot Hedlin and Constable! . 



176 Il990l ). Aj are damping factors, Yid and R^ are weighting matrices to give dif- 

177 ferent elements of data misfit and model roughness vectors equal weights, 

178 Arrij is the change in model parameters for the i-th iteration and Adj is the 

179 data misfit vector containing the difference between calculated and observed 

180 apparent resistivities. Since the Ad values may extend over several orders of 

181 magnitude, logarithmic differences are employed. 

182 Equation H] is solved iteratively until either the root-mean square (RMS) of 

183 the data misfit vector Adj does not change significantly after an inversion 

184 step and/or it becomes smaller than the measurement accuracy. The weight- 

185 ing matrices R^ and Rm are predefined, and default values were chosen for 

186 A,- . 
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1B7 2.2 Ensembles 



IBS Inversion problems for geoelectrical surveys are usually ill-posed, mixed deter- 

189 mined problems. If the errors in data acquisition and in the inversion process 

190 would be known quantitatively, the optimum model and its error distribution 

191 could be determined exactly. Measurement errors often can only be estimated, 

192 and further discrepancies may be introduced during inversion, especially if an 

193 inversion code is used that does not rigorously optimise for a given error esti- 

194 mate. Additionally, inverted models can be plagued by possibly large inversion 

195 artifacts depending e.g. on resistivity contrasts. 

196 2.2.1 Building Ensembles 

197 Consequently, it might not be sufficient to analyse only the optimum model 

198 (i.e. the model with the smallest data misfit), but to compute a range of 

199 possible models addressing the inherent variability of the inversion process. 

200 By randomly varying the inversion parameter set and creating an ensemble 

201 of possible inversion models, the whole parameter space and thus the possible 

202 model range is explored. 

203 For the RES2DINV code used here, the selected parameters are listed in Table 

204 [U The table also includes for each parameter the range from which a value was 

205 automatically and randomly selected. The parameter selection encompasses 

206 the use of smoothness constrained and robust inversions as well as two mixed 

207 formulations with a robust constraint applied only on the data, and one with 

208 a robust constraint applied only on the model. Further variations address 

209 the regularisation, e.g. the damping factor, where an initial damping factor 

210 Xstart and the maximum damping factor Xmax are varied. For most variations. 
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211 the maximum damping factor Xmax is kept at Xmax = 10 • Xstart- Additional 

212 variations include the reduction of side block effects, the ratio of vertical to 

213 horizontal smoothness filtering and the use of the first iteration step model 

214 as a reference model for the further iterations instead of using the average of 

215 resistivities. 

216 It should be noted that this choice of variations is specific for the software used 

217 in this study. However, the idea can easily be transferred to similar inversion 

218 approaches. 

219 Almost all inversions resulted in inverted models with RMS errors smaller 

220 than 4% as can be expected from adding 3% artifical noise to the data set. 

221 Some single inversions, however, resulted in a larger RMS error. In section [3l 

222 both, inversion models with RMS < 4% and > 4%, will be included to keep 

223 the ensembles balanced. 



224 2.2.2 Clustering 



225 Each ensemble is created as a set of 50 different inversion models and then 



19881). 50 



226 regrouped using a fc-means clustering algorithm (iDubes and Jainl . 

227 models have been chosen arbitrarily as a compromise between computational 

228 efficiency and the necessity to generate a sufficiently large ensemble for clus- 

229 tering. The fc-means clustering method starts with a collection of genes; here 

230 a gene is a row of all block resistivities of one model. The distance d between 

231 two genes is calculated as a Pearson correlation 



CTr 



Vi-y 



(5) 



233 where x is the average of values in gene x and is the standard deviation 



234 of these values ( lEisen et al. 



19981 ). The fc- means clustering starts with a user 
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235 decision on the number of clusters to be created, then randomly assigns each 

236 gene to a cluster. For each cluster, the average model is created, then each gene 

237 is assigned to the cluster it has the smallest distance from. These last steps 

238 are repeated until an optimal solution is found. At least two r uns creating 



239 the s ame optimal solution are needed to reach a reliable solution (lEisen et al. 



240 Il998l ). In this study, we used five clusters to generate a sufficiently large cluster 

241 variability while ensuring that the number of ensemble members per cluster 

242 is not too small. 



243 2.3 Applicability 

244 All five steps presented here form an analysis cycle for a synthetic case study 

245 investigating the reliability of resistivity quantificaton for shallow subsurface 

246 water processes. For application to field cases, it is possible to create and 

247 analyse a simplified synthetic representation of the actual site following the 

248 five steps above o rapplay only the ensemble and clustering steps t determine 

249 the spread of possible inversion results. 



250 3 Synthetic case studies 



251 All test cases studied here are based on a simple two layer mediur n represent- 



252 i ng th e structure of a full-scale dike model described in detail by 



Rings et al 



253 ( 120081 ). Although synthetic data sets are employed to distinctively focus on 

254 specific anomalies, the material parameters were obtained from real observa- 

255 tion. Hydraulic parameters, following the van Genuchte n-Mualem par a meter - 



256 isation, were determined in a laboratory experiment by 



Scheuermann! ([20051) 
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257 For obtaining the parameters of the overburden, we used an inversion proce- 

258 dure supported by the HYDRUS software. As no direct measurements of water 

259 content in the over 



Durden were available, a rainfall experiment described in 
Scheuermann! (120051 ) was simulated. Pressure head measurements in the sand, 
but directly below the overburden, were used to invert the hydraulic param- 
eters of the overburden. The resulting parameters are listed in Table [2] for 
the two different materials. Meteorological data from the permanent station 
Karlsruhe-Nordwest (German y) were used as forcing . The Penman formula. 



265 calibrated for grass cover by 



Doorenbos and Pruitt 



(Il977l ). was applied to 



these data to retrieve values for potential evapotranspiration. Combined with 
measured precipitation rates, these values have been used as daily averages 
for simulations of 210 days based on measurements in 2001. 
Based on this two layer medium, two generic cases representing different ide- 
alised case studies were created: The first case study simulates a defective 
sealing, where an infiltration plume of water is generated in the sand layer. In 
the second case study, a rectangular, hydraulically resistive anomaly is placed 
in the sand underneath the organic overburden. 



274 3.1 First Case: Defective Sealing 

275 The first case is based on the idea of a crack in a dike sealing. Damaged 

276 sealings are critical, as even through small cracks, large amounts of water can 

277 infiltrate. 

278 In this hypothetical case, water infiltrates through an otherwise sealed off 

279 surface through one crack. The sealing is considered to be invisible to the 

280 geoelectrical survey. 
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2B1 3.1.1 Water Simulation 



2B2 In HYDRUS, the sealing is modeled as a no-flow boundary, and the crack 

2B3 has an atmospheric boundary and is filled with sand material. The simula- 

284 tion results show water infiltrating through the crack into the sand where it 

285 diffuses into a sinking plume. The water content does not change outside of 

286 the plume (Fig. [2]). Three characteristic states of the simulated results can 

287 be identified: dry state (Fig. [2^), infiltration state (the plume begins to form 

288 in the sand, Fig. Eb) and the diffusion state (Fig. [2b), where the center of 

289 the plume has propagated into the sand and the top layer is already drying. 

290 The transfer from water content to resistivities was done by assuming a dry 

291 state resistivity of p = 400 Qm for the overburden and p = 5000 Qm for the 

292 sand and applying Eq. [3] with satur ation exponent n = 2 for the overbur- 



293 den and n = 1.164 for the sand (see 



Rings et al. 



20081 ). During infiltration 



294 and diffusion, this results in a minimal resistivity in the plume of p = 2000 Qm. 



296 3.1.2 Forward- Inverse Cycle 

297 Figure [3] shows three standard (robust) inversion models for the three states 

298 of water percolation. A complete Wenner-Schlumberger array with electrode 

299 separation 0.5 m has been simulated in the forward modeling. In the dry state, 

300 the crack is clearly visible. In the infiltration state, the infiltrating plume is 

301 characterised through a distinct lower resistivity than the background, while 

302 in the diffusion state the inversion did not sufficiently resolve the shape of the 

303 plume. 

304 To analyse the dependence of the inversion results on the resistivity con- 
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305 trast between the plume and the host material, the plume resistivity was 

306 increased or lowered in steps of 250 Qm around the minimal plume resistivity 

307 of 2000 Qm. A total of nine models with plume resistivity ranging from 1000 

308 to 3000 flm were explored, while the background resistivity stayed constant 

309 at 5000 Qm. 

310 Generally, the resistivity of an anomaly p~nom 

Panom = min{pi} (6) 



312 for all model blocks i below the overburden. The misfit in the anomaly's 

313 resistivity Apm is the difference between the resistivity of the anomaly in the 

314 generic p'nom^gen and inverted model p~nom,inv- 

Ap (7) 



316 For this case study, p^nom corresponds to the resistivity in the center of the 

317 plume. Figure H] shows the results of the forward-inverse cycle as Apm vs the 

318 resistivity contrast. While the error in resistivity quantification is smallest 

319 for the orginal contrast of 4:10, smaller and higher contrasts both result in 

320 increasingly larger Apm- 

321 Apm is slightly smaller in the infiltration state. In the diffusion state, the 

322 center of the plume has sunk to greater depth, where the reduced sensitivity 

323 of ERT may be the reason for a less accurate quantification. 



324 3.1.3 Ensemble 

325 The inversion ensemble for the case of the defective sealing and the diffusion 

326 state is shown in Figure [5l All models within the ensemble detected the over- 
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327 burden with the damaged seahng, but the model parts below this overburden 

328 show different features. In the first cluster, fi-sloped artifacts appear to the 

329 side of the plume with equal resistivity as the plume itself. In the second 

330 cluster the artifacts appear as well, but have comparably higher resistivity, so 

331 that the plume appears as a distinct feature. In the third cluster, both plume 

332 and fi-sloped artifacts are roughly in the same resistivity range, but have 

333 a higher resistivity than in cluster 1. The fourth cluster comprises strongly 

334 damped models where the plume is mostly visible. The last cluster shows 

335 models where the plume is clearly visible, with comparably better contrast, 

336 but mostly the vertical extent of the plume feature is overestimated. 

337 To comprehend the ensemble results in a simple way, averaged models of each 

338 cluster are shown in Figure O As the clustering process already involves av- 

339 eraging, this is a valid method. In Figure [6l the mean models for each of 

340 the clusters of the ensemble shown in Figure [5] are now listed according to 

341 the number of cluster members. It must be noted that the smallest cluster 

342 contains only 3 models, whereas the largest cluster contains almost half the 

343 models of the ensemble. The average RMS error of each cluster is below 4%. 

344 The most prominent feature retrieved in all models is the two-layered struc- 

345 ture, which can be observed in all five clusters. This structure is present even 

346 in clusters where the damping is strong enough to nearly hide the plume 

347 anomaly. When comparing clusters 3-5 to the strongly damped inversion re- 

348 suits in cluster 1, the typical fi-sloped structure can be identified as an artifact 

349 at the lateral boundaries of the plume. Compared to the standard model (0), 

350 the cluster averages allow a much better identification of features, even though 

351 some interpretational experience or a priori knowledge is needed to distinguish 

352 between real anomalies (cluster 5) and artifacts (cluster 4). 
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353 3.2 Second Case: Hydmulically Resistive Anomaly 

354 In the second case, the accuracy of resistivity quantification for a rectangular, 

355 hydrauhcally resistive anomaly placed below the organic overburden is stud- 

356 ied. First, a soil model with an organic overburden and an anomaly at 0.55 

357 m depth was created in HYDRUS. To represent the hydrauhcally resistive 

358 material of the anomaly, the same material as for the organic overburden was 

359 used. Then, multiple versions of this model were created with slightly differ- 

360 ent geometries. Table [3] shows the differences between the respective models, 

361 which will be explained in the following. 

362 

363 3.2.1 Water Simulation 

364 In the simulation of water movement, a dry state, an infiltration state and a 

365 diffusion state were identified as characteristic states of an infiltrating water 

366 front. In the dry state (Fig. [7^), the soil is completely free of water. In the 

367 infiltration state (Fig.lTb), the water front is propagating into the volume. The 

368 hydrauhcally resistive anomaly causes water to impound on top, only slowly 

369 infiltrating into the anomaly. In the diffusion state (Fig. [Tt), the infiltration 

370 front has reached the bottom boundary of the model, and the organic overbur- 

371 den and parts of the sand directly below are beginning to dry. The anomaly 

372 is filled with water that infiltrates into the sand beneath. 

373 Analysis of the quality of water content estimation through ERT was con- 

374 ducted for a variety of models and electrode configurations based on the three 

375 states of water percolation in Figure [71 To study the influence of contrasting 

376 resistivities at the surface, models with and without an organic overburden 
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377 were used for simulation. In addition, the depth of the anomaly was varied in 

378 steps of 0.2 m with the upper boundary at 0.35 m to 1.15 m depth. To examine 

379 the effect of electrode configuration, two different electrode arrays (complete 
3B0 Wenner-Schlumberger and Dipole-Dipole arrays) with an electrode spacing of 
381 0.5 m were used for each model (Table [3]). 

382 

383 3.2.2 Forward-Inverse Cycle 

384 Inspection of the inverted models (Fig. [HI right column) shows that the rectan- 

385 gular shape of the anomaly cannot be exactly retrieved. Determination of an 

386 average resistivity of the anomaly would be dependent on an arbitrary deter- 

387 mination of anomaly borders. It is also not possible to determine the average 

388 resistivity at the actual position of the anomaly, since the perceived depth of 

389 the anomaly is greater than the actual depth. 

390 In the following, results for the different models shown in Table [3] will be 

391 compared regarding Ap^ (Eq- E|), which now corresponds to the (minimal) re- 

392 sistivity of the anomaly. Figure [3 shows Ap^ as a function of anomaly depth. 

393 For the Wenner-Schlumberger array, Apm increases with anomaly depth, reach- 

394 ing up to 2-3 times the expected value. A much better estimate is obtained if 

395 no organic overburden is present (gray curves). For these cases, better quan- 

396 tifications of p~nom possible and Ap^ increases only slightly with depth. In 

397 the diffusion state, significantly smaller errors occur compared to other states 

398 of water percolation, especially in the presence of an organic overburden. 

399 As can be seen in Figure El the error in depth resolution is rather large. If 

400 an organic overburden is present, the thickness of this layer is overestimated, 

401 causing a shift in the vertical position of the anomaly of 0.3 to 0.4 m. It was 
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402 also observed that at greater depths, the position stays approximately the 

403 same for an anomaly expected at 0.75 m to 1.15 m depth. Again, in the case 

404 of a model without an organic overburden, the higher sensitivity due to higher 

405 resistivities near the surface makes better depth determination possible. 

406 For models simulated with the Dipole-Dipole array, errors for models with or- 

407 ganic overburden are significantly smaller than for the Wenner-Schlumberger 

408 array. However, the Dipole-Dipole array was shown to be very sensitive to 

409 noise and disturbances at the surface (like a stone pathway), to a point were 

410 measurements taken using this array could not be interpreted with the avail- 

411 able inversion routines. 

412 As a measure of the quality of the inversion, a simple criterion containing the 

413 model misfit M as the sum of all errors has been applied: 



415 where Fi is the i-ih. model block of the inversion domain discretisation. 

416 Comparison of M for the different states of water percolation (Fig. [TOj) shows 

417 that the diffusion state gives significantly better results. In this state, the misfit 

418 below the organic overburden and to the sides of the anomaly is much smaller, 

419 additionally the depth of the anomaly and overburden are better resolved. 

420 Figure [TT] shows the spatial error distribution for each state. In the dry state, 

421 the biggest errors stem from an overestimated thickness of the overburden, 

422 which also entails further mispositioning of the anomaly. The anomaly itself is 

423 also vertically elongated, leading to considerable errors in the lower parts. In 

424 the diffusion state (Fig. [11] (b)), the resistivity contrast between overburden 

425 and wet sand is much smaller, due to a) the sand having a reduced resistivity 

426 as it is more saturated with water and b) the overburden being dryer as in 
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427 the previous states, resulting in a higher resistivity. As a consequence of this 

428 reduced resistivity contrast, the errors resulting from an incorrect overburden 

429 thickness are reduced as well. 

430 

431 3.2.3 Ensemble 

432 For the case of the hydraulically resistive anomaly, the random set of param- 

433 eters is applied to generic models of all three different states of water perco- 

434 lation. To assure comparability, the random parameter set stays the same for 

435 each of the three models. 

436 A model with an anomaly at 0.75 m depth was used, including an organic 

437 overburden and using the Wenner-Schlumberger array. For each state, an en- 

438 semble of 50 inverted models was created. For simplification, only mean cluster 

439 members are shown. Figure fT2] shows the five clusters per ensemble with the re- 

440 spective number of ensemble members. The respective Ap^ is listed in Table HI 

441 

442 • Dry State: The rectangular shape of the anomaly is retrieved variably well, 

443 but for the models 4 and 5, where the thickness of the anomaly is smaller, 

444 a strong overestimation of resistivities is present in the lower part of the 

445 model (> 7000 Vim instead of 5000 Vtm). The resistivity of the anomaly 
Panom ^uch too high for all five models. For models 1 and 2 that contain 

447 most of the ensemble members, the anomaly is vertically elongated. 

448 • Infiltration State: In four models, the shape of the anomaly has been re- 

449 trieved quite well, but for model 10, two zones of minimal resistivity have 

450 been detected rather than the rectangular shape. In all models, Ap^ is very 
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451 large compared to the expected resistivity of the anomaly p = 65 Qm. Again, 

452 models 9 and 10 (same inversion parameters as model 4 and 5) overestimate 

453 the background resistivity at greater depth. 

454 • Diffusion State: The resistivity of the anomaly is detected with lower resis- 

455 tivity as in the infiltration state, closer to the expected resistivity of 45 Qm. 

456 Again, in model 13 and 15, the strong inversion artifact is present near the 

457 bottom coinciding with the shape of the anomaly being retrieved quite well. 

458 These artifacts are not present in model 11, 12 and 14, where the anomaly is 

459 vertically elongated. Model 15 presents a mixed case of a slightly elongated 

460 anomaly and an artifact of smaller extent than in model 13. 

461 TableHlshows, sorted for the cluster representatives, the misfits in the anomaly's 

462 resistivity. While it is apparent that the errors are large in each case, they are 

463 again considerably smaller for the diffusion state. 

464 4 Discussion and Conclusion 

465 The ability of electrical resistivity tomography to accurately determine resis- 

466 tivity distributions was examined. A two-step model approach was used to 

467 create synthetic data sets. It comprises the modeling of soil water movement 

468 for synthetic soil data sets and a transfer into a model of generic resistivi- 

469 ties using a petrophysical relation. A forward-inverse cycle is used evaluate 

470 how well the geophysical inversion scheme can reconstruct the given soil data 

471 set and its water content. An ensemble and clustering approach is proposed 

472 because a single model deduced as the optimal model does not necessarily 

473 reproduce the expected resistivities accurately. 

474 The methods were applied to two case studies of simple soil models based on a 
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475 two-layered structure reproducing field observations. The first case simulates 

476 the infiltration of water through a cracked surfical sealing, and the second a 

477 hydraulically resistive anomaly in a sand layer. 

478 Key results of the forward-inverse modeling in this study include: 

479 • In the presence of large resistivity contrasts, e.g. a conductive organic over- 

480 burden, the retrieval of accurate resistivity values beneath this layer using 

481 the regularisation based inversion method applied in this study is not possi- 

482 ble. However, if the volume is monitored at various stages of water percola- 

483 tion, the retrieval quality can differ. Especially in the diffusion state, much 

484 better accuracy was possible. 

485 • The model misfit increases with depth, as the sensitivity of the inversion 

486 model to the data decreases. 

487 • In the absence of an organic overburden, a much better quantification is 

488 possible because of a lower resistivity contrast. 

489 • The numerical study showed that a Dipolc-Dipolc array provides more accu- 

490 rate inversions than the Wenner-Schlumberger array. However, in practical 

491 applications, it has to be ensured that the signal-to-noise ratio is sufficiently 

492 large. 

493 As a consequence, an ensemble approach was introduced that creates multiple 

494 inversion models for one data set by randomly choosing the inversion parame- 

495 ters from the possible (and numerically plausible) parameter space. By using 

496 clustering methods, averaged models representing different clusters in the en- 

497 semble can be created and compared. Key results of the ensemble approach 

498 include: 

499 



22 



500 • Clustering of ensemble members allows an evaluation of the different possi- 

501 ble models that fit the data. Areas likely to be plagued by artifacts can be 

502 identified and the reliability of standard inverted models can be evaluated. 

503 • However, the quantification of resistivities is not considerably improved by 

504 ensembles. For example, it became apparent that resistivities retrieved with 

505 smaller misfits in one region can coincide with larger artifacts in other re- 

506 gions. 

507 • The clustering of ensembles allows an overview of the ensemble, without 

508 losing information about the ensemble. 

509 The ideas of the approaches presented here can easily be adapted to different 

510 models and inversion methods. For the specific inversion process with regular- 

511 isation used in this study, it can be concluded that a reliable quantification of 

512 resistivity values is not possible. The use of additional information, e.g. within 

513 a framework aiming at directly inverting or calculating hydrological proper- 

514 ties from collected data sets that not only contain resistivity measurements, 

515 but also data about the flow conditions, e.g. meteorological data, should be 

516 considered. 
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Constraint on the data 


robust or smooth 


Constraint on the model 


robust or smooth 


Initial damping Aj 


0.01 to 1 


Minimal damping 


0.05Ai to 0.2Ai 


Convergence limit 


1% to 9% 


Maximal number of iterations 


3, 5 or 15 


Vertical to horizontal regularisation 


0.25 to 4 


Increase of damping with depth 


1.0 to 2.0 


Reduce effect of 
side blocks 


none, slight, 
severe, very severe 


Higher damping for first layer 


yes or no 



Table 1 

Parameter space of inversion parameters used for ensemble calculations. 

629 rithms, covergence analysis and numerical comparisons. Statistical Com- 

630 putations 9, 907-921. 

631 Zhou, B., Dahlin, T., 2003. Properties and effects of measurement errors on 

632 2d resistivity imaging surveying. Near Surface Geophysics 1, 105-117. 
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Material 




ds 


a 


n 


Ks [m/d] 


Sand 


0.045 


0.361 


4 


2.2 


17.28 


Overburden 


0.067 


0.45 


5.23 


2.67 


0.225 



Table 2 

Soil parameters for the van Genuchten-Mualem parameterisation. 9r is the residual 
water content, 9s is the volumetric water content at full saturation, a and n are 
parameters connected to the pore radii, and Ks is the hydraulic conductivity at 
saturation. 



Organic 
overburden 


with 


without 








Depth of 
the anomaly [m] 


0.35 


0.55 


0.75 


0.95 


1.15 


State of percolation 


dry 


infiltration 


diffusion 






ERT array 


WS 


DD 









Table 3 

Parameter variation for different soil models, infiltration state and measurement 
geometry 
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Cluster DRY 


Cluster INFILTRATION 


Cluster DIFFUSION 


1) 381 


6) 156 


11) 49 


2) 613 


7) 294 


12) 41 


3) 228 


8) 314 


13) 86 


4) 1068 


9) 73 


14) 62 


5) 608 


10) 292 


15) 228 



Table 4 

Misfit for the cluster representative shown in Figure [12] (misfits in Om) . 
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Water 
Movement 
in Soil 





Fig. 1. Charts visualizing the methodological steps involved in this study. Above: 
Steps in the forward-inverse cycle. Below: Steps in the ensemble method. 
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Fig. 2. Defective sealing, characteristic states of water percolation, (a) Dry State 
(b) Infiltration State (c) Diffusion State. 
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Fig. 3. Inverted models for the case of the defective seahng. (a) Dry State (b) 
Infiltration State (c) Diffusion State. 
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Fig. 4. Misfit of tlie anomaly for the case of the defective seahng, shown is the 
resistivity contrast as the ratio plume divided by host material (x-axis) versus Ap^ 
of the the inverted model (y-axis). 
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Fig. 5. Clustered ensemble with 50 possible models for the diffusion state of the case 
of the defective sealing with an infiltration plume. The domains of the 5 clusters 
are indicated by numbers and dividing lines. 
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Fig. 6. Standard model (0) and averaged cluster models (1-5) for the case of the 
defective sealing. In contrast to Figure [5l the clusters are sorted in descending order 
by the number of ensemble members. 
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Fig. 7. States of the simulation of water movement through a model with a hy- 
draulically resistive anomaly (rectangular block marked with thick black outline). 
The layer boundary between organic overburden and sand is marked with a thin 
horizontal line, (a) Dry State (b) Infiltration State (c) Diffusion State. 
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Fig. 8. Generic and inverted models for the anomaly (Wenner-Schlumberger ar- 
ray), (a) Dry State (b) Infiltration State (c) Diffusion State. The black rectangle 
in the right column marks the location of the anomaly in the left column, the thin 
horizontal line marks the layer boundary between organic overburden and sand. 




Fig. 9. Apra for cases with organic overburden (black lines) and without (gray lines). 
Top row: Survey with Wenner-Schlumberger, bottom row: with Dipole-Dipole. The 
left column shows the dry state, the middle column the infiltration state and the 
right column the diffusion state. 
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Fig. 10. Cumulative block misfits for the three stages of water percolation. Shown is 
the logarithm of the sum of all errors M for Wenner-Schlumberger and Dipole-Dipole 
arrays and with or without an organic overburden. 
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Fig. 12. Averaged cluster representatives for the resistive anomaly in the dry (left), 
infiltration (middle) and diffusion (right) state. 
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